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ABSTRACT 


The ability to predict failure loads in notched composite laminates is a requirement in a variety of 
structural design circumstances. A complicating factor is the development of a zone of damaged 
material around the notch tip. The objective of this study was to develop a computational 
technique that simulates progressive damage growth around a notch in a manner that allows the 
prediction of failure over a wide range of notch sizes. This was accomplished through the use of 
a relatively simple, nonlocal damage model that incorporates strain-softening. This model was 
implemented in a two-dimensional finite element program. Calculations were performed for two 
different laminates with various notch sizes under tensile loading, and the calculations were found 
to correlate well with experimental results. 
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1 INTRODUCTION 


The design of structural components made of composite materials is heavily influenced by damage 
tolerance requirements. The problem of predicting failure in notched laminates has been the 
subject of numerous studies. A review of those occurring before 1985 can be found in Awerbuch 
and Madhukar. 1 All failure theories that have been developed contain empirical parameters such 
as a characteristic length. Unfortunately, these parameters are not always true material properties 
because they appear to be functions of notch size; i.e., the parameter value that leads to a correct 
failure prediction for small notches does not usually work well with large notches. 2 

In a composite material a zone of damage of considerable influence is known to develop in 
advance of a crack. Recent research 3 indicates that the damage growth in the vicinity of the crack 
tip manifests itself in the form of strain-softening of the material. Strain-softening models have 
received a great deal of attention for describing the fracture of concrete and other materials where 
microstructure has a strong influence on macroscopic properties. In general, a macroscopic 
description of damage (eg., distributed cracking) is reflected in a constitutive model that exhibits 
a decrease of stress with increasing strain beyond some critical strain e d as shown in Fig. 1 . The 
incorporation of strain-softening into a finite element analysis based on classical plasticity theory 
results in calculations that are mesh sensitive. This occurs because as the mesh is refined, there is 
a tendency for the damage zone to localize to a zero volume. This leads to the prediction of 
structural failure with zero energy dissipation, which is physically impossible. Numerous 
techniques have been proposed to address this issue. 

Two of these techniques that hold promise for modeling damage growth in composite 
structures are the discrete crack model and the nonlocal damage model. The discrete crack model 
avoids the zero energy dissipation problem by using a stress-displacement law over the damage 
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zone rather than a stress-strain law. These models were pioneered by Hillerborg et al . 4 who used 
a Dugdale-Barrenblat 5,6 type model to represent the damage zone ahead of the crack. Backlund 
and Aronsson 7,8 extended this concept to composite materials. This approach is straightforward 
but lacks generality because the crack path must be assumed beforehand. Alternatively, the 
nonlocal damage model offers the same advantages of the discrete crack model but avoids its 
shortcomings. 

Rudimentary damage mechanics based upon local constitutive relations can be used to 
successfully model strain softening. Better yet, damage evolution based upon the concept of a 
nonlocal continuum prevents unacceptable localization of the damage. The basic idea in a 
nonlocal continuum is that the stress at a point in a body is dependent on the strains of neighbor- 
ing points in addition to the strain at the point. In its most general form, this theory, which can be 
tied to a statistical analysis of heterogeneous materials, can lead to nearly intractable mathematical 
complexity. Bazant and Chang 9 modified this theory and incorporated it into a finite element 
analysis that could handle strain-softening without spurious mesh sensitivity. Unfortunately, this 
approach proved to be inconvenient in practical applications because it required an imbricated 
system of elements overlapping one another. Recently, Bazant and Pijaudier-Cabot 10 removed 
this difficulty with the introduction of a nonlocal damage theory which incorporates nonlocal 
principles into continuum damage mechanics. In this theory only the damage parameter is 
considered to be nonlocal (i.e., a function of the strain averaged from a certain neighborhood of a 
point) while all other variables are local. Bazant 11 justified the validity of this approach using a 
simplified micromechanical analysis to show that fracturing strain due to damage is the result of 
the release of stored energy from a finite size microcrack neighborhood. Bazant and Lin 12 
incorporated the nonlocal damage concept into a two-dimensional finite element analysis and 
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applied it successfully in predicting failure of concrete. DeVree et al. 13 also developed a nonlocal 
damage finite element analysis for isotropic materials. 

The objective of this study was to develop a technique that simulates progressive damage 
growth around a notch in a composite laminate and that can be used to predict failure. This was 
accomplished by applying a relatively simple nonlocal damage model to fiber-reinforced compos- 
ites. This model was implemented in a two-dimensional finite element program. Calculations 
were performed, and the results were compared to data from fracture tests. 

2 STRESS-STRAIN RELATIONS FOR A DAMAGED LAMINATE 

The analysis of a composite laminate will not be done on a ply-by-ply basis, but rather it will 
simply be considered as a homogeneous material with orthotropic material properties. Treating 
the material as being in a state of plane stress, the stress-strain relations in principal material 
coordinates (xj and x 2 ) for undamaged material can be written as 

{€} = [C]- 1 {0} (1) 

where 

{e} T = [e, e 2 Yn] ( 2 ) 

{o} T =[o, o 2 t i2 ] (3) 

C„ = Ej/(1 - v 12 v 21 ) , C 22 = E^l - v I2 v 21 ) , C 33 = G 12 

C]2 — c 21 — v 12 E 2 /( 1 - v 12 v 21 ) , Cj 3 = C^j — C 31 = C 32 — 0 . (4) 

When damage occurs (microcracking, etc.) the effective load-carrying area of the material is 
reduced. We introduce the concept of an effective stress 14 d = o/(l-D) to account for this area 
reduction. The quantity D is the damage variable which ranges from 0 (no damage) to 1 
(development of a macrocrack). For simplicity we assume that damage develops independently in 
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the Xj and x 2 directions. For this orthotropic damage development, the effective stress can be 
expressed in general form as 15 

{5} = [M] {o} (5) 


where 


{d} T = [d, d 2 t I2 ] 


M n = 1/(1 -Di) , M ^ = 1/(1 -D 2 ) , M 33 = 1/fi^ f^D 2 

M 12 — M|j M23 ^^21 ~ ^31 = M32 = 0 • 


( 6 ) 


(7) 


We now take the same approach as Chow and Wang 15 and use the elastic energy equivalence 
concept that states that the complementary elastic energy for a damaged material is in the same 
form as that of an undamaged material except that the stress is replaced by the effective stress in 
the energy formulation, i.e., 


W = - {o} T [C ]' 1 {0} 
2 


( 8 ) 


= | {o) T (M] T [C]-‘ [M] { 0 } 


( 9 ) 


The stress-strain equation for a damaged material can be written as 
{ £ I = = [M] T [C] » [M] {0} 


or 


3(0} 


(£( = [C]' 1 {O} 


where 


[C] > = [M] T [C] 1 [M] 
Alternatively, we can write 


( 10 ) 


( 11 ) 


( 12 ) 


5 



<°> - [C] {£} 


(13) 


where 


C„ = Ej (1-D,) 2 / (l-v 12 v 21 ) 

C J2 = C 21 = VjjEj (1-D,) (1-D 2 ) / (l-v 12 v 21 ) 

C 22 = E 2 (1-D 2 ) 2 / (l-v, 2 v 21 ) (14) 

C 33 = Gj 2 (1-D,) (1-D 2 ) 

= E 23 = C 31 = C 32 = 0 . 

In the finite element formulation, we may have the material coordinates x, and x 2 at some 


angle 0 relative to the global coordinates x and y. We can transform the stresses and strains 


through the usual transformation relations to get 
{o} = [ TJ {o'} 
and 

{e} = [T c ] {€'} 

where 


{o'} T = 

\Px °y ^xy] 



{e'} T = 

[Cx Yxy] 




COS 2 0 

sin 2 0 

2sin0 cos0 

[TJ = 

sin 2 0 

cos 2 0 

-2sin0 cos0 


-sin0 cos0 

sin0 cos0 

cos 2 0 - sin 2 0 


cos 2 0 

sin 2 0 

sin0 cos0 

[TJ = 

sin 2 0 

cos 2 0 

-sin0 cos0 


-2sin0 cos0 

2sin0 cos0 cos 2 -sin 2 0 


(15) 

(16) 

(17) 

(18) 

(19) 

( 20 ) 
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Combining eqns (13), (15), and (16) gives 


{O'} = [C'J {€'} 


( 21 ) 


where 


[C'j = [TJ 1 [C] [TJ . 


( 22 ) 


3 DAMAGE DEVELOPMENT 


In the nonlocal damage model, the damage parameters D t and D 2 are assumed to be 


functions of the nonlocal strains e x and e 2 defined as 


e,(x,y) = 


e 2 (x,y) ■ 


— 7— : ff a (^ _x > tf - y) e 1 (^» T l) d£ dr| 
A r (x,y) JJ a 

— — - ff a(^-x,Ti-y)e 2 (^ti)d^dil 

A r (x,y) JJ a 


(23) 

(24) 


where 

A r (x,y) = ff a(£-x, r\ -y) d£ d r\ (25) 

and a(£-x, rj-y) is a weight function which we have chosen in the form suggested by Bazant, 16 i.e., 

a(£-x,n-y) = ji - 


(^-x) 2 + (n -y) 2 
0.8256 { 2 


(26) 


where { is the characteristic length for the material, and a equals zero if the quantity within 
brackets becomes negative. The relationship between the damage parameter and the nonlocal 
strain is determined from the stress-strain curve. Consider the case of a uniaxial load o 0 in the x t 
direction. For this case eqn (1 1) gives 
e, = o„/E,(l -D,f . 
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Since e 2 is uniform, the nonlocal strain e, is simply equal to e, for all points except those in a 
boundary layer around the edges which will be ignored. Solving eqn (27) for D t gives 

D, = 1 - (<J 0 /E,I,)'' 2 . (28) 

The shape of the stress-strain curve after the initiation of damage (i.e., the relationship between o L 
and e, for > e ld ) determines how Dj is related to Cj. This ultimately leads to an explicit 
mathematical representation for Dj as 

Di = fj (ij) (29) 

A similar relationship can be found for D 2 , i.e., 

D 2 = f 2 (i 2 ). (30) 

Equations (21), (23), (24), (29), and (30) completely describe the stress-strain behavior of the 
laminate. 

4 FINITE ELEMENT FORMULATION 

To develop a finite element formulation for progressive damage analysis, we begin with 
the principle of virtual work 17 

/■ («'} T (o'|dV = f (u}*{f}ds (31) 

J V J s' 

where {e'} is the strain associated with the virtual displacement {u}, {u} s is the virtual displace- 
ment of the surface of the body, and {f} is the traction on the surface of the body. We will 
develop this analysis for an 8-node quadrilateral element. Using the usual shape functions of this 
element, 1 * we can express the displacement [u] m within an element m in terms of the nodal 
displacement {U} as 

{u} m = [L] m {U} . (32) 
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Applying the two-dimensional strain-displacement relations to eqn (32), we arrive at the following 


relation 18 between strain and nodal displacement 

( e }m = [B]m {U} . (33) 

Substituting eqns (32) and (33) into eqn (31) gives 

E / P»£<o').dV. - E f (34) 

m V m m J S m 

Substituting eqn (21) into eqn (34) gives 

(E f [B]I[C'l„[B]„dvJ {U } = E / Ms.{0„dS„ (35) 

V « JV n> ) m J S ro 

We now set 

M = 'LL [B£[C']„[B].dV„, (36) 

m Jv m 

W =E f ILfi.{0.dS.. (37) 

m Js m 

Thus 

[K] (U} = {R} (38) 


where [K] is the stiffness matrix and {R} is the generalized nodal load matrix. After damage 
initiates, eqn (38) represents a nonlinear system of algebraic equations because [K] is a function of 
D, and D 2 which depend on the nonlocal strains. The nonlocal strains in eqns (23) and (24) are 
evaluated numerically using Gaussian integration. Equation (38) can be solved iteratively using 
the Newton-Raphson method. 17 

Based on the above analysis, a computer program was developed and calculations were 
performed for the case of a center-cracked plate under tension, as shown in Fig. 2. It was found 
that convergence difficulties arose for materials whose stress-strain curves had softening regimes 
with steep slopes. To overcome these difficulties we made use of a procedure based on the 
concept of viscous relaxation. 19 We introduced a small amount of viscous damping into the 
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analysis and treated the problem dynamically (without inertia effects) so that eqn (38) was 
replaced by 

[C d ] {U} + [K] {U} = {R} (39) 

where [C d ] is the damping matrix and {U} is the time derivative of {U}. To solve the above 
differential equations numerically, we used the trapezoidal rule of time integration, 17 i.e., 

HAt {U} = ‘{U} + At ( ‘{U} + t+At {U} )/2 (40) 

where the superscript in front of the variable indicates the time at which it is evaluated. For the 
damping matrix we employed Rayleigh damping, i.e., 

[CJ = P [K] . (41) 

We found that a value of P = 0.01 sec' 1 for a time step At = 1 sec eliminated the conver- 
gence problems and also reproduced results that were within a fraction of a percent of previous 
static calculations. 

5 RESULTS 

Calculations were performed for a 13 -ply laminate of graphite/epoxy with a 
[45/-45/90/0/60/-60/90/-60/60/0/90/-45/45] lay-up. The laminate stiffness properties were E x 
=36.3 GPa, Ey =61.5 GPa, = 0.28, and G^ =20.7 Gpa.. The softening portions of the stress- 
strain curves were represented by exponential functions: 

°X = (42) 

°y - E y € <iy e ' VV ‘*' > (43) 

with = 0.01 101, a* = 2000, = 0.00984, ay = 600, and a characteristic length ? = 0.762 cm. 

The stress-strain curve for the y-direction is shown in Fig. 1. The values for and e^, came 
from tensile tests. 2 Measuring the softening properties in a standard tensile test is extremely 
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difficult to do. 20 Therefore, these values were chosen in an inverse manner to match data from 
fracture tests described next. 

A center-cracked plate under tension as shown in Fig. 2 was analyzed with the program 
for a crack length 2a = 6.35 cm. and a specimen width W = 25.4 cm. An applied load o 0 was 
simulated by gradually displacing the top boundary while holding the bottom boundary fixed as 
would be done in a tensile test. Load versus displacement was monitored during the calculation, 
and when the load reached its peak value and began to decline, the specimen was assumed to fail. 
Calculations were performed for three different quarter-symmetry meshes as shown in Fig. 3 
where MESH 1 is a coarse mesh, MESH 2 is a moderately fine mesh, and MESH 3 is a very fine 
mesh. A reflecting algorithm had to be used to correctly calculate nonlocal strains in elements 
near a symmetry boundary. The load versus displacement results for these three meshes are 
shown in Fig. 4. It can be seen that the solution converges with mesh refinement, and the 
difference in failure load between MESH 2 and MESH 3 is less than 0.5 percent. The calculations 
were repeated for various crack lengths between 0 and 30.5 cm. with W/2a = 4. A plot of 
nominal failure stress (applied force divided by total cross-sectional area) versus crack length is 
shown in Fig. 5. Also shown are experimental results from fracture tests. 2 It can be seen that the 
theory is capable of representing the failure load over a wide range of crack lengths. 

The extent of the damage zone (i.e., material that has passed into the strain-softening 
regime) was also monitored during the calculation. As expected, the damage zone originated near 
the crack tip and grew primarily in the horizontal direction away from the tip. The distance from 
the crack tip to the outer edge of the damage zone as a function of load for three different crack 
lengths is shown in Fig. 6. In each case damage does not initiate until the applied load is very 
close to the failure load. After damage initiates, it grows rapidly. 
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The calculations just described were repeated for a 15-ply laminate of graphite/epoxy with 
a [-45/45/0/90/-3 0/30/-75/0/75/30/-30/90/0/45/-45] lay-up. The laminate stiffness properties 
were E x =56.0 Gpa, Ey=46.5 Gpa, v^O.342, and 0^=19. 6 Gpa. Again, the softening portions of 
the stress-strain curves were represented by exponential functions with €^=0.0072, a x =20, 

€ d y = 0 0072, ay=20, and a characteristic length {=0.254 cm. The stress-strain curve for the y- 
direction for this laminate is shown in Fig. 7. Comparing this to Fig. 1, we observe that the 
softening portion of the curve for this laminate is much less steep than that of the 13-ply laminate. 

A plot of nominal failure stress versus crack length for the 15-ply laminate is shown in Fig. 
8 along with the corresponding experimental results. Again, the theory has performed well in 
representing the failure load over a wide range of crack lengths. The distance to the outer edge of 
the damage zone as a function of load for three crack lengths is shown in Fig. 9. Comparing this 
to the results in Fig. 6, we observe that the growth of the damage zone with applied load for this 
laminate is considerably more gradual than that for the 13-ply laminate. Overall, the 15-ply 
laminate exhibits less brittle behavior than does the 13 -ply laminate as one would expect from the 
softening portion of their stress-strain curves. 

6 CONCLUSION 

In this study a progressive damage model for predicting failure loads in notched composite 
laminates was developed. This was based on relatively simple, nonlocal damage mechanics 
incorporating strain-softening. A two-dimensional finite element program was developed, and 
calculations were performed for two different laminates under tensile loading. One laminate 
exhibited “brittle” strain-softening response, and the other exhibited “ductile” strain-softening 
response. The theory was shown to accurately predict failure of these laminates over a wide 
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range of notch sizes. This approach clearly shows promise as a design tool for assessing damage 
tolerance because it has the potential for studying the response of a structure after a crack has 
begun to propagate. However, a number of issues need further study. These include the choice 
of mathematical function used to represent the strain-softening portion of the stress-strain curve 
and the effect of flexure. These topics will be addressed in future research. 
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